#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gmp.h>
#include <time.h>

int main(int argc, char *argv[]) {
  /* Váriaveis para contar o tempo */
  clock_t begin, end;

  /* Váriaveis do algoritmo */
  mpf_t an, an1;
  mpf_t bn, bn1;
  mpf_t tn, tn1;
  mpf_t pn, pn1;
  mpf_t pi;

  FILE *pi_file = fopen("pi_gauss_sequencial.txt","w+");

  if(pi_file == NULL) {
    return EXIT_FAILURE;
  }

  long N; /* Número de iterações */
  long i;

  if(argc != 2) {
    printf("\nUsage: %s [n_iter].\n",argv[0]);
    exit(-1);
  } else {
    sscanf(argv[1],"%ld",&N);
  }

  mpf_set_default_prec(40000000); /* Precisão necessária para calcular 10 milhões de dígitos */

  begin = clock();
  /* Inicialização das váriaveis */
  mpf_init_set_ui(an,1);

  mpf_init(bn);
  mpf_sqrt_ui(bn,2);
  mpf_ui_div(bn,1,bn);

  mpf_init_set_str(tn,"0.25",10);

  mpf_init_set_ui(pn,1);

  mpf_init(an1);
  mpf_init(bn1);
  mpf_init(tn1);
  mpf_init(pn1);

  mpf_init(pi);

  printf("Inciando cálculo do PI pelo método de Gauss-Legendre com %ld iterações.\n",N);


  for(i=0;i<N;i++) {   
   
    mpf_init(an1);
    //an1 = (an + bn)/2.0;
    mpf_add(an1,an,bn);
    mpf_div_ui(an1,an1,2);

    mpf_init(bn1);
    //bn1 = pow(an*bn,1.0/2.0);
    mpf_mul(bn1,an,bn);
    mpf_sqrt(bn1,bn1);

    mpf_init(tn1);
    //tn1 = tn - pn*(an - an1)*(an - an1);
    mpf_sub(tn1,an,an1);  // tn1 = an - an1
    mpf_mul(tn1,tn1,tn1); // tn1 = (an - an1)^2
    mpf_mul(tn1,tn1,pn);  // tn1 = pn*(an - an1)^2    
    mpf_sub(tn1,tn,tn1);  

    mpf_init(pn1);
    //pn1 = 2*pn;
    mpf_mul_ui(pn1,pn,2);

    //an = an1;
    mpf_set(an,an1);
    //bn = bn1;
    mpf_set(bn,bn1);
    //tn = tn1;
    mpf_set(tn,tn1);
    //pn = pn1;
    mpf_set(pn,pn1);

    printf("Iteração: %ld\n",i);
  
    mpf_add(pi,an,bn);   // pi = an + bn
    mpf_pow_ui(pi,pi,2); // pi = (an + bn)^2
    mpf_div_ui(pi,pi,4); // pi = (an + bn)^2 / 4
    mpf_div(pi,pi,tn);   // pi = (an + bn)^2 / 4 / tn*/
    
    mpf_out_str(pi_file,10,0,pi);
    fprintf(pi_file,"\n");    
  }


  
  end = clock();
  printf("\nTempo do cálculo: %lf segundos\n",(double)((end - begin)/CLOCKS_PER_SEC));

  fclose(pi_file);
  mpf_clears(pi,an,bn,tn,pn,an1,bn1,tn1,pn1,NULL);

}
